
import pandas as pd
import matplotlib.pyplot as plt #if using matplotlib
import plotly.express as px #if using plotly
import geopandas as gpd
import warnings

import momepy

import warnings

import geopandas as gpd
import libpysal
import momepy
import osmnx as ox
import pandas as pd
import os
import pandas as pd
import geopandas as gpd
import geopy
from geopy.geocoders import Nominatim
from geopy.extra.rate_limiter import RateLimiter
import matplotlib.pyplot as plt

import tqdm
from tqdm._tqdm_notebook import tqdm_notebook


area = pd.read_csv('areafinal.csv')


area = area[['MSOA11CD', 'MSOA11NM',
       'MSOA21CD', 'MSOA21NM', 'LAD22CD', 'LAD22NM', 'CTRY21CD', 'CTRY21NM',
       'TCITY15CD', 'TCITY15NM']]


area.drop_duplicates(subset=['MSOA11CD'], keep='first')


area.head()


MSOA = gpd.read_file('MSOA.geojson')


MSOA.plot()

<AxesSubplot:>


MSOA = pd.merge(MSOA,area,how='left',left_on='MSOA11CD',right_on='MSOA11CD')


MSOA.tail()


MSOA['centroid'] = MSOA.centroid

C:\Users\neetm\AppData\Local\Temp\ipykernel_11612\2953065240.py:1: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  MSOA['centroid'] = MSOA.centroid


print(MSOA['centroid'][0])

POINT (-0.0923889468524288 51.514381323272744)


MSOA = MSOA.set_geometry('geometry')


MSOA["area"] = MSOA.area

C:\Users\neetm\AppData\Local\Temp\ipykernel_11612\1553680121.py:1: UserWarning: Geometry is in a geographic CRS. Results from 'area' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  MSOA["area"] = MSOA.area


MSOA = MSOA[['OBJECTID', 'MSOA11CD', 'MSOA11NM_x','MSOA21CD', 'MSOA21NM', 'LAD22CD', 'LAD22NM',
       'CTRY21CD', 'CTRY21NM', 'TCITY15CD', 'TCITY15NM','BNG_E', 'BNG_N', 'LONG_', 'LAT',
       'Shape_Leng', 'Shape__Area', 'Shape__Length', 'geometry','centroid', 'area']]


MSOA = MSOA.rename({'MSOA11NM_x':'MSOA11NM'})


MSOA = MSOA[MSOA['MSOA11CD'].str[0]=='E']


MSOA.plot("area", legend=True)

<AxesSubplot:>


MSOA.explore("area", legend=False)


eng_df = MSOA


place = 'London, United Kingdom'
local_crs = 27700


import geopandas as gpd
gpd.tools.geocode(place).explore()


x=0
MSOA0 = MSOA.iloc[x].geometry


buildings = ox.geometries.geometries_from_polygon(MSOA.iloc[x].geometry, tags={'building':True})
buildings = buildings[buildings.geom_type == 'Polygon'].reset_index()
buildings = buildings[["geometry"]].to_crs(local_crs)
buildings["uID"] = range(len(buildings))
buildings.head()


blg_area = momepy.Area(buildings)
buildings['area'] = blg_area.series
buildings.head()


buildings.explore()


f, ax = plt.subplots(figsize=(10, 10))
buildings.plot('area', ax=ax, legend=True, scheme='quantiles', cmap='Blues',
               legend_kwds={'loc': 'lower left'})
ax.set_axis_off()
plt.show()


f, ax = plt.subplots(figsize=(10, 10))
buildings.plot(ax=ax)
ax.set_axis_off()
plt.show()


blg_ERI = momepy.EquivalentRectangularIndex(buildings)


blg_ERI = momepy.EquivalentRectangularIndex(buildings, areas='area')
buildings['eri'] = blg_ERI.series


f, ax = plt.subplots(figsize=(10, 10))
buildings.plot('eri', ax=ax, legend=True)
ax.set_axis_off()
plt.show()


buildings.head()


osm_graph = ox.graph_from_polygon(MSOA.iloc[x].geometry, network_type='drive')
osm_graph = ox.projection.project_graph(osm_graph, to_crs=local_crs)
streets = ox.graph_to_gdfs(
    osm_graph,
    nodes=False,
    edges=True,
    node_geometry=False,
    fill_edge_geometry=True)
streets = momepy.remove_false_nodes(streets)
streets = streets[["geometry"]]
streets["nID"] = range(len(streets))
streets.head()

D:\Anaconda\envs\geo_env\lib\site-packages\geopandas\array.py:1406: UserWarning: CRS not set for some of the concatenation inputs. Setting output's CRS as OSGB36 / British National Grid (the single non-null crs provided).
  warnings.warn(


streets.explore()


limit = momepy.buffered_limit(buildings, 100)
tessellation = momepy.Tessellation(buildings, "uID", limit, verbose=False, segment=1)
tessellation = tessellation.tessellation

D:\Anaconda\envs\geo_env\lib\site-packages\momepy\elements.py:384: UserWarning: Tessellation does not fully match buildings. 1 element(s) collapsed during generation - unique_id: {1236}
  warnings.warn(
D:\Anaconda\envs\geo_env\lib\site-packages\momepy\elements.py:395: UserWarning: Tessellation contains MultiPolygon elements. Initial objects should be edited. unique_id of affected elements: [1799, 405, 655, 1719, 446, 654, 663, 445, 317, 638, 633, 634, 640, 631, 796, 799, 815, 551, 671, 1786, 723, 552, 174, 1796, 1459, 1787, 1797, 1782, 842, 1450, 621, 1785, 1403, 291, 188, 1829, 1830, 181, 1461, 189, 947, 1466, 179, 1744, 629, 844, 1922, 1924, 851, 857, 860, 343, 1165, 1736, 644, 1628, 274, 1278, 1938, 591, 1937, 1276, 1270, 134, 1267, 1269, 142, 1629, 1811, 270, 1886, 1810, 1123, 1122, 1129, 1130, 1124, 1112, 1232, 1930, 229, 1912, 1929, 1923, 1458, 1317, 1045, 87, 297, 1365, 1368, 1360, 95, 1500, 570, 573, 571, 1058, 98, 1869, 1792]
  warnings.warn(


tessellation.explore()


tessellation.plot()

<AxesSubplot:>


tes_area = momepy.Area(tessellation)
tessellation['area'] = tes_area.series


tessellation.head()


tessellation.head()


f, ax = plt.subplots(figsize=(10, 10))
tessellation.plot(ax=ax, column='area', legend=True, scheme='quantiles', k=10, cmap='viridis')
buildings.plot(ax=ax, color='white', alpha=0.5)
ax.set_axis_off()
plt.show()


tes_cwa = momepy.CompactnessWeightedAxis(tessellation)
tessellation['cwa'] = tes_cwa.series


f, ax = plt.subplots(figsize=(10, 10))
tessellation.plot(ax=ax, column='cwa', legend=True, scheme='quantiles', k=10, cmap='Greens_r')
ax.set_axis_off()
plt.show()


buildings = buildings.sjoin_nearest(streets, max_distance=1000, how="left")
buildings.head()


buildings = buildings.drop_duplicates("uID").drop(columns="index_right")
tessellation = tessellation.merge(buildings[['uID', 'nID']], on='uID', how='left')


buildings["area"] = buildings.area
tessellation["area"] = tessellation.area
streets["length"] = streets.length
buildings['eri'] = momepy.EquivalentRectangularIndex(buildings).series
buildings['elongation'] = momepy.Elongation(buildings).series
tessellation['convexity'] = momepy.Convexity(tessellation).series
streets["linearity"] = momepy.Linearity(streets).series
buildings['perimeter'] = momepy.Perimeter(buildings).series
buildings['wall_length'] = momepy.PerimeterWall(buildings).series
streets['length_neighbours'] = momepy.SegmentsLength(streets, mean=True).mean

Calculating spatial weights...
Spatial weights ready...

  0%|          | 0/1947 [00:00<?, ?it/s]

  0%|          | 0/1947 [00:00<?, ?it/s]

Calculating spatial weights...
Spatial weights ready...

  0%|          | 0/1221 [00:00<?, ?it/s]


sw = libpysal.weights.DistanceBand.from_dataframe(tessellation,
                                                  threshold=100,
                                                  silence_warnings=True,
                                                  ids='uID')
tessellation['mean_area'] = momepy.AverageCharacter(tessellation,
                                                    values='area',
                                                    spatial_weights=sw,
                                                    unique_id='uID').mean

  0%|          | 0/1946 [00:00<?, ?it/s]


fig, ax = plt.subplots(1, 2, figsize=(24, 12))

buildings.plot("eri", ax=ax[0], scheme="natural_breaks", legend=True)
buildings.plot("elongation", ax=ax[1], scheme="natural_breaks", legend=True)

ax[0].set_axis_off()
ax[1].set_axis_off()

D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=8.
  warnings.warn(
D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=8.
  warnings.warn(


fig, ax = plt.subplots(1, 2, figsize=(24, 12))

tessellation.plot("convexity", ax=ax[0], scheme="natural_breaks", legend=True)
streets.plot("linearity", ax=ax[1], scheme="natural_breaks", legend=True)

ax[0].set_axis_off()
ax[1].set_axis_off()

D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=8.
  warnings.warn(
D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=5.
  warnings.warn(


buildings["shared_walls"] = momepy.SharedWallsRatio(buildings).series
buildings.plot("shared_walls", figsize=(12, 12), scheme="natural_breaks", legend=True).set_axis_off()

D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=8.
  warnings.warn(


queen_1 = libpysal.weights.contiguity.Queen.from_dataframe(tessellation, ids="uID", silence_warnings=True)


tessellation["neighbors"] = momepy.Neighbors(tessellation, queen_1, "uID", weighted=True, verbose=False).series
tessellation["covered_area"] = momepy.CoveredArea(tessellation, queen_1, "uID", verbose=False).series

with warnings.catch_warnings():
    warnings.simplefilter("ignore")

    buildings["neighbor_distance"] = momepy.NeighborDistance(buildings, queen_1, "uID", verbose=False).series


fig, ax = plt.subplots(1, 2, figsize=(24, 12))

buildings.plot("neighbor_distance", ax=ax[0], scheme="natural_breaks", legend=True)
tessellation.plot("covered_area", ax=ax[1], scheme="natural_breaks", legend=True)

ax[0].set_axis_off()
ax[1].set_axis_off()

D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=8.
  warnings.warn(
D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=8.
  warnings.warn(


queen_3 = momepy.sw_high(k=3, weights=queen_1)
buildings_q1 = libpysal.weights.contiguity.Queen.from_dataframe(buildings, silence_warnings=True)

buildings['interbuilding_distance'] = momepy.MeanInterbuildingDistance(buildings, queen_1, 'uID', queen_3, verbose=False).series
buildings['adjacency'] = momepy.BuildingAdjacency(buildings, queen_3, 'uID', buildings_q1, verbose=False).series


fig, ax = plt.subplots(1, 2, figsize=(24, 12))

buildings.plot("interbuilding_distance", ax=ax[0], scheme="natural_breaks", legend=True)
buildings.plot("adjacency", ax=ax[1], scheme="natural_breaks", legend=True)

ax[0].set_axis_off()
ax[1].set_axis_off()

D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=8.
  warnings.warn(
D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=8.
  warnings.warn(


profile = momepy.StreetProfile(streets, buildings)
streets["width"] = profile.w
streets["width_deviation"] = profile.wd
streets["openness"] = profile.o


fig, ax = plt.subplots(1, 3, figsize=(24, 12))

streets.plot("width", ax=ax[0], scheme="natural_breaks", legend=True)
streets.plot("width_deviation", ax=ax[1], scheme="natural_breaks", legend=True)
streets.plot("openness", ax=ax[2], scheme="natural_breaks", legend=True)

ax[0].set_axis_off()
ax[1].set_axis_off()
ax[2].set_axis_off()

D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=5.
  warnings.warn(
D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=5.
  warnings.warn(
D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=5.
  warnings.warn(


buildings.head()


tessellation['car'] = momepy.AreaRatio(tessellation, buildings, 'area', 'area', 'uID').series
tessellation.plot("car", figsize=(12, 12), vmin=0, vmax=1, legend=True).set_axis_off()


graph = momepy.gdf_to_nx(streets)
graph = momepy.node_degree(graph)
graph = momepy.closeness_centrality(graph, radius=400, distance="mm_len")
graph = momepy.meshedness(graph, radius=400, distance="mm_len")
nodes, streets = momepy.nx_to_gdf(graph)

  0%|          | 0/657 [00:00<?, ?it/s]

  0%|          | 0/657 [00:00<?, ?it/s]


fig, ax = plt.subplots(1, 3, figsize=(24, 12))

nodes.plot("degree", ax=ax[0], scheme="natural_breaks", legend=True, markersize=1)
nodes.plot("closeness", ax=ax[1], scheme="natural_breaks", legend=True, markersize=1, legend_kwds={"fmt": "{:.6f}"})
nodes.plot("meshedness", ax=ax[2], scheme="natural_breaks", legend=True, markersize=1)

ax[0].set_axis_off()
ax[1].set_axis_off()
ax[2].set_axis_off()

D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=3.
  warnings.warn(
D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=3.
  warnings.warn(
D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=3.
  warnings.warn(


buildings["nodeID"] = momepy.get_node_id(buildings, nodes, streets, "nodeID", "nID")

  0%|          | 0/1947 [00:00<?, ?it/s]


tessellation.head()


merged = tessellation.merge(buildings.drop(columns=['nID', 'geometry']), on='uID')
merged = merged.merge(streets.drop(columns='geometry'), on='nID', how='left')
merged = merged.merge(nodes.drop(columns='geometry'), on='nodeID', how='left')


merged.columns

Index(['uID', 'geometry', 'area_x', 'nID', 'convexity', 'mean_area',
       'neighbors', 'covered_area', 'car', 'area_y', 'eri', 'elongation',
       'perimeter', 'wall_length', 'shared_walls', 'neighbor_distance',
       'interbuilding_distance', 'adjacency', 'nodeID', 'length', 'linearity',
       'length_neighbours', 'width', 'width_deviation', 'openness', 'mm_len',
       'node_start', 'node_end', 'degree', 'closeness', 'meshedness'],
      dtype='object')


blocks = momepy.Blocks(tessellation, streets, buildings, 'bID', 'uID')
blocks.blocks.head()


tessellation.head()


blocks.blocks.explore()


limit = momepy.buffered_limit(buildings)


limit


sw = momepy.sw_high(k=3, gdf=tessellation, ids='uID')
tessellation['covered'] = momepy.CoveredArea(tessellation, sw, 'uID').series

  0%|          | 0/1946 [00:00<?, ?it/s]


tessellation.explore()


y = 462
print(MSOA.iloc[y].MSOA11CD)

E02000481


import osmnx


buildings = ox.geometries.geometries_from_polygon(MSOA.iloc[y].geometry, tags={'building':True})
buildings = buildings[buildings.geom_type == 'Polygon'].reset_index()
buildings = buildings[["geometry"]].to_crs(local_crs)
buildings["uID"] = range(len(buildings))
buildings.head()


buildings.explore()


buildings['MSOA11CD'] = MSOA.iloc[y].MSOA11CD
osm_graph = osmnx.graph_from_polygon(MSOA.iloc[y].geometry, network_type='drive')
osm_graph = osmnx.projection.project_graph(osm_graph, to_crs= local_crs)
streets = osmnx.graph_to_gdfs(osm_graph,nodes=False,edges=True,node_geometry=False,fill_edge_geometry=True)


streets.explore()


fig, ax = plt.subplots(1, 2, figsize=(24, 12))

buildings.plot("eri", ax=ax[0], scheme="natural_breaks", legend=True)
buildings.plot("elongation", ax=ax[1], scheme="natural_breaks", legend=True)

ax[0].set_axis_off()
ax[1].set_axis_off()

D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=5.
  warnings.warn(
D:\Anaconda\envs\geo_env\lib\site-packages\sklearn\cluster\_kmeans.py:1334: UserWarning: KMeans is known to have a memory leak on Windows with MKL, when there are less chunks than available threads. You can avoid it by setting the environment variable OMP_NUM_THREADS=5.
  warnings.warn(


streets = momepy.remove_false_nodes(streets)
streets = streets[["geometry"]]
streets["nID"] = range(len(streets))
    
limit = momepy.buffered_limit(buildings, 100)
tessellation = momepy.Tessellation(buildings, "uID", limit, verbose=False, segment=1)
tessellation = tessellation.tessellation
buildings = buildings.sjoin_nearest(streets, max_distance=1000, how="left")
buildings = buildings.drop_duplicates("uID").drop(columns="index_right")
tessellation = tessellation.merge(buildings[['uID', 'nID']], on='uID', how='left')
buildings["area"] = buildings.area
tessellation["area"] = tessellation.area
streets["length"] = streets.length
    
buildings['eri'] = momepy.EquivalentRectangularIndex(buildings).series
buildings['elongation'] = momepy.Elongation(buildings).series
tessellation['convexity'] = momepy.Convexity(tessellation).series
streets["linearity"] = momepy.Linearity(streets).series
    
buildings["shared_walls"] = momepy.SharedWallsRatio(buildings).series
queen_1 = libpysal.weights.contiguity.Queen.from_dataframe(tessellation, ids="uID", silence_warnings=True)
tessellation["neighbors"] = momepy.Neighbors(tessellation, queen_1, "uID", weighted=True, verbose=False).series
tessellation["covered_area"] = momepy.CoveredArea(tessellation, queen_1, "uID", verbose=False).series
    
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    buildings["neighbor_distance"] = momepy.NeighborDistance(buildings, queen_1, "uID", verbose=False).series
    queen_3 = momepy.sw_high(k=3, weights=queen_1)
    buildings_q1 = libpysal.weights.contiguity.Queen.from_dataframe(buildings, silence_warnings=True)
    buildings['interbuilding_distance'] = momepy.MeanInterbuildingDistance(buildings, queen_1, 'uID', queen_3, verbose=False).series
    buildings['adjacency'] = momepy.BuildingAdjacency(buildings, queen_3, 'uID', buildings_q1, verbose=False).series
    profile = momepy.StreetProfile(streets, buildings)
    streets["width"] = profile.w
    streets["width_deviation"] = profile.wd
    streets["openness"] = profile.o
    tessellation['car'] = momepy.AreaRatio(tessellation, buildings, 'area', 'area', 'uID').series
    graph = momepy.gdf_to_nx(streets)
    graph = momepy.node_degree(graph)
    graph = momepy.closeness_centrality(graph, radius=400, distance="mm_len")
    graph = momepy.meshedness(graph, radius=400, distance="mm_len")
    nodes, streets = momepy.nx_to_gdf(graph)
    buildings["nodeID"] = momepy.get_node_id(buildings, nodes, streets, "nodeID", "nID")
    merged = tessellation.merge(buildings.drop(columns=['nID', 'geometry']), on='uID')
    merged = merged.merge(streets.drop(columns='geometry'), on='nID', how='left')
    merged = merged.merge(nodes.drop(columns='geometry'), on='nodeID', how='left')

D:\Anaconda\envs\geo_env\lib\site-packages\geopandas\array.py:1406: UserWarning: CRS not set for some of the concatenation inputs. Setting output's CRS as OSGB36 / British National Grid (the single non-null crs provided).
  warnings.warn(

  0%|          | 0/116 [00:00<?, ?it/s]

  0%|          | 0/116 [00:00<?, ?it/s]

  0%|          | 0/55 [00:00<?, ?it/s]


merged = merged.drop(['uID','geometry','nID'],axis=1)


y = pd.DataFrame(merged.mean()).transpose()
display(y)

C:\Users\neetm\AppData\Local\Temp\ipykernel_11612\2010261618.py:1: FutureWarning: Dropping of nuisance columns in DataFrame reductions (with 'numeric_only=None') is deprecated; in a future version this will raise TypeError.  Select only valid columns before calling the reduction.
  y = pd.DataFrame(merged.mean()).transpose()


x2 = 737
MSOA.iloc[x2].geometry


buildings = ox.geometries.geometries_from_polygon(MSOA.iloc[x2].geometry, tags={'building':True})
buildings = buildings[buildings.geom_type == 'Polygon'].reset_index()
buildings = buildings[["geometry"]].to_crs(local_crs)
buildings["uID"] = range(len(buildings))
buildings.head()


buildings.explore()


buildings['MSOA11CD'] = MSOA.iloc[x2].MSOA11CD
osm_graph = osmnx.graph_from_polygon(MSOA.iloc[x2].geometry, network_type='drive')
osm_graph = osmnx.projection.project_graph(osm_graph, to_crs= local_crs)
streets = osmnx.graph_to_gdfs(osm_graph,nodes=False,edges=True,node_geometry=False,fill_edge_geometry=True)


streets.explore()


streets = momepy.remove_false_nodes(streets)
streets = streets[["geometry"]]
streets["nID"] = range(len(streets))
    
limit = momepy.buffered_limit(buildings, 100)
tessellation = momepy.Tessellation(buildings, "uID", limit, verbose=False, segment=1)
tessellation = tessellation.tessellation
buildings = buildings.sjoin_nearest(streets, max_distance=1000, how="left")
buildings = buildings.drop_duplicates("uID").drop(columns="index_right")
tessellation = tessellation.merge(buildings[['uID', 'nID']], on='uID', how='left')
buildings["area"] = buildings.area
tessellation["area"] = tessellation.area
streets["length"] = streets.length
    
buildings['eri'] = momepy.EquivalentRectangularIndex(buildings).series
buildings['elongation'] = momepy.Elongation(buildings).series
tessellation['convexity'] = momepy.Convexity(tessellation).series
streets["linearity"] = momepy.Linearity(streets).series
    
buildings["shared_walls"] = momepy.SharedWallsRatio(buildings).series
queen_1 = libpysal.weights.contiguity.Queen.from_dataframe(tessellation, ids="uID", silence_warnings=True)
tessellation["neighbors"] = momepy.Neighbors(tessellation, queen_1, "uID", weighted=True, verbose=False).series
tessellation["covered_area"] = momepy.CoveredArea(tessellation, queen_1, "uID", verbose=False).series
    
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    buildings["neighbor_distance"] = momepy.NeighborDistance(buildings, queen_1, "uID", verbose=False).series
    queen_3 = momepy.sw_high(k=3, weights=queen_1)
    buildings_q1 = libpysal.weights.contiguity.Queen.from_dataframe(buildings, silence_warnings=True)
    buildings['interbuilding_distance'] = momepy.MeanInterbuildingDistance(buildings, queen_1, 'uID', queen_3, verbose=False).series
    buildings['adjacency'] = momepy.BuildingAdjacency(buildings, queen_3, 'uID', buildings_q1, verbose=False).series
    profile = momepy.StreetProfile(streets, buildings)
    streets["width"] = profile.w
    streets["width_deviation"] = profile.wd
    streets["openness"] = profile.o
    tessellation['car'] = momepy.AreaRatio(tessellation, buildings, 'area', 'area', 'uID').series
    graph = momepy.gdf_to_nx(streets)
    graph = momepy.node_degree(graph)
    graph = momepy.closeness_centrality(graph, radius=400, distance="mm_len")
    graph = momepy.meshedness(graph, radius=400, distance="mm_len")
    nodes, streets = momepy.nx_to_gdf(graph)
    buildings["nodeID"] = momepy.get_node_id(buildings, nodes, streets, "nodeID", "nID")
    merged = tessellation.merge(buildings.drop(columns=['nID', 'geometry']), on='uID')
    merged = merged.merge(streets.drop(columns='geometry'), on='nID', how='left')
    merged = merged.merge(nodes.drop(columns='geometry'), on='nodeID', how='left')

D:\Anaconda\envs\geo_env\lib\site-packages\geopandas\array.py:1406: UserWarning: CRS not set for some of the concatenation inputs. Setting output's CRS as OSGB36 / British National Grid (the single non-null crs provided).
  warnings.warn(

  0%|          | 0/88 [00:00<?, ?it/s]

  0%|          | 0/88 [00:00<?, ?it/s]

  0%|          | 0/36 [00:00<?, ?it/s]


merged = merged.drop(['uID','geometry','nID'],axis=1)
x2 = pd.DataFrame(merged.mean()).transpose()
display(x2)

C:\Users\neetm\AppData\Local\Temp\ipykernel_11612\4099531547.py:2: FutureWarning: Dropping of nuisance columns in DataFrame reductions (with 'numeric_only=None') is deprecated; in a future version this will raise TypeError.  Select only valid columns before calling the reduction.
  x2 = pd.DataFrame(merged.mean()).transpose()


last = 850
MSOA.iloc[last].geometry


buildings = ox.geometries.geometries_from_polygon(MSOA.iloc[last].geometry, tags={'building':True})
buildings = buildings[buildings.geom_type == 'Polygon'].reset_index()
buildings = buildings[["geometry"]].to_crs(local_crs)
buildings["uID"] = range(len(buildings))
buildings.head()


osmnx.projection.project_geometry(MSOA.iloc[last].geometry)

(<shapely.geometry.polygon.Polygon at 0x179ce327a30>,
 <Derived Projected CRS: +proj=utm +zone=30 +ellps=WGS84 +datum=WGS84 +unit ...>
 Name: unknown
 Axis Info [cartesian]:
 - E[east]: Easting (metre)
 - N[north]: Northing (metre)
 Area of Use:
 - undefined
 Coordinate Operation:
 - name: UTM zone 30N
 - method: Transverse Mercator
 Datum: World Geodetic System 1984
 - Ellipsoid: WGS 84
 - Prime Meridian: Greenwich)


buildings.explore()


buildings['MSOA11CD'] = MSOA.iloc[last].MSOA11CD
osm_graph = osmnx.graph_from_polygon(MSOA.iloc[last].geometry, network_type='drive')
osm_graph = osmnx.projection.project_graph(osm_graph, to_crs= local_crs)
streets = osmnx.graph_to_gdfs(osm_graph,nodes=False,edges=True,node_geometry=False,fill_edge_geometry=True)


streets.explore()


streets = momepy.remove_false_nodes(streets)
streets = streets[["geometry"]]
streets["nID"] = range(len(streets))
    
limit = momepy.buffered_limit(buildings, 100)
tessellation = momepy.Tessellation(buildings, "uID", limit, verbose=False, segment=1)
tessellation = tessellation.tessellation
buildings = buildings.sjoin_nearest(streets, max_distance=1000, how="left")
buildings = buildings.drop_duplicates("uID").drop(columns="index_right")
tessellation = tessellation.merge(buildings[['uID', 'nID']], on='uID', how='left')
buildings["area"] = buildings.area
tessellation["area"] = tessellation.area
streets["length"] = streets.length
    
buildings['eri'] = momepy.EquivalentRectangularIndex(buildings).series
buildings['elongation'] = momepy.Elongation(buildings).series
tessellation['convexity'] = momepy.Convexity(tessellation).series
streets["linearity"] = momepy.Linearity(streets).series
    
buildings["shared_walls"] = momepy.SharedWallsRatio(buildings).series
queen_1 = libpysal.weights.contiguity.Queen.from_dataframe(tessellation, ids="uID", silence_warnings=True)
tessellation["neighbors"] = momepy.Neighbors(tessellation, queen_1, "uID", weighted=True, verbose=False).series
tessellation["covered_area"] = momepy.CoveredArea(tessellation, queen_1, "uID", verbose=False).series
    
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    buildings["neighbor_distance"] = momepy.NeighborDistance(buildings, queen_1, "uID", verbose=False).series
    queen_3 = momepy.sw_high(k=3, weights=queen_1)
    buildings_q1 = libpysal.weights.contiguity.Queen.from_dataframe(buildings, silence_warnings=True)
    buildings['interbuilding_distance'] = momepy.MeanInterbuildingDistance(buildings, queen_1, 'uID', queen_3, verbose=False).series
    buildings['adjacency'] = momepy.BuildingAdjacency(buildings, queen_3, 'uID', buildings_q1, verbose=False).series
    profile = momepy.StreetProfile(streets, buildings)
    streets["width"] = profile.w
    streets["width_deviation"] = profile.wd
    streets["openness"] = profile.o
    tessellation['car'] = momepy.AreaRatio(tessellation, buildings, 'area', 'area', 'uID').series
    graph = momepy.gdf_to_nx(streets)
    graph = momepy.node_degree(graph)
    graph = momepy.closeness_centrality(graph, radius=400, distance="mm_len")
    graph = momepy.meshedness(graph, radius=400, distance="mm_len")
    nodes, streets = momepy.nx_to_gdf(graph)
    buildings["nodeID"] = momepy.get_node_id(buildings, nodes, streets, "nodeID", "nID")
    merged = tessellation.merge(buildings.drop(columns=['nID', 'geometry']), on='uID')
    merged = merged.merge(streets.drop(columns='geometry'), on='nID', how='left')
    merged = merged.merge(nodes.drop(columns='geometry'), on='nodeID', how='left')

D:\Anaconda\envs\geo_env\lib\site-packages\geopandas\array.py:1406: UserWarning: CRS not set for some of the concatenation inputs. Setting output's CRS as OSGB36 / British National Grid (the single non-null crs provided).
  warnings.warn(

  0%|          | 0/131 [00:00<?, ?it/s]

  0%|          | 0/131 [00:00<?, ?it/s]

  0%|          | 0/1075 [00:00<?, ?it/s]


tessellation.explore()


import pandas as pd
area = pd.read_excel('area_final.xlsx')
area = area.drop('Unnamed: 0',axis=1)
csv = pd.read_csv('msoacity.csv')
csv.info()
area.head()
area.columns
area = area[['MSOA11CD', 'MSOA11NM', 'LSOA21CD', 'LSOA21NM',
       'MSOA21CD', 'MSOA21NM', 'LAD22CD', 'LAD22NM', 'CTRY21CD', 'CTRY21NM',
       'TCITY15CD', 'TCITY15NM']]
area = area.drop_duplicates(subset=['MSOA11CD'], keep='first')

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 6791 entries, 0 to 6790
Data columns (total 5 columns):
 #   Column     Non-Null Count  Dtype 
---  ------     --------------  ----- 
 0   MSOA11CD   6791 non-null   object
 1   MSOA11NM   6791 non-null   object
 2   TCITY15CD  3526 non-null   object
 3   TCITY15NM  4196 non-null   object
 4   FID        6791 non-null   int64 
dtypes: int64(1), object(4)
memory usage: 265.4+ KB


area.columns

Index(['MSOA11CD', 'MSOA11NM', 'LSOA21CD', 'LSOA21NM', 'MSOA21CD', 'MSOA21NM',
       'LAD22CD', 'LAD22NM', 'CTRY21CD', 'CTRY21NM', 'TCITY15CD', 'TCITY15NM'],
      dtype='object')

	MSOA11CD	MSOA11NM	MSOA21CD	MSOA21NM	LAD22CD	LAD22NM	CTRY21CD	CTRY21NM	TCITY15CD	TCITY15NM
0	E02000001	City of London 001	E02000001	City of London 001	E09000001	City of London	E92000001	England	J01000055	London
1	E02000002	Barking and Dagenham 001	E02000002	Barking and Dagenham 001	E09000002	Barking and Dagenham	E92000001	England	J01000055	London
2	E02000003	Barking and Dagenham 002	E02000003	Barking and Dagenham 002	E09000002	Barking and Dagenham	E92000001	England	J01000055	London
3	E02000004	Barking and Dagenham 003	E02000004	Barking and Dagenham 003	E09000002	Barking and Dagenham	E92000001	England	J01000055	London
4	E02000005	Barking and Dagenham 004	E02000005	Barking and Dagenham 004	E09000002	Barking and Dagenham	E92000001	England	J01000055	London
...	...	...	...	...	...	...	...	...	...	...
6786	E02006930	Greenwich 037	E02006930	Greenwich 037	E09000011	Greenwich	E92000001	England	J01000055	London
6787	E02006931	Greenwich 038	E02006931	Greenwich 038	E09000011	Greenwich	E92000001	England	J01000055	London
6788	E02006932	Liverpool 060	E02006932	Liverpool 060	E08000012	Liverpool	E92000001	England	J01000054	Liverpool
6789	E02006933	Liverpool 061	E02006933	Liverpool 061	E08000012	Liverpool	E92000001	England	J01000054	Liverpool
6790	E02006934	Liverpool 062	E02006934	Liverpool 062	E08000012	Liverpool	E92000001	England	J01000054	Liverpool

	MSOA11CD	MSOA11NM	MSOA21CD	MSOA21NM	LAD22CD	LAD22NM	CTRY21CD	CTRY21NM	TCITY15CD	TCITY15NM
0	E02000001	City of London 001	E02000001	City of London 001	E09000001	City of London	E92000001	England	J01000055	London
1	E02000002	Barking and Dagenham 001	E02000002	Barking and Dagenham 001	E09000002	Barking and Dagenham	E92000001	England	J01000055	London
2	E02000003	Barking and Dagenham 002	E02000003	Barking and Dagenham 002	E09000002	Barking and Dagenham	E92000001	England	J01000055	London
3	E02000004	Barking and Dagenham 003	E02000004	Barking and Dagenham 003	E09000002	Barking and Dagenham	E92000001	England	J01000055	London
4	E02000005	Barking and Dagenham 004	E02000005	Barking and Dagenham 004	E09000002	Barking and Dagenham	E92000001	England	J01000055	London

	OBJECTID	MSOA11CD	MSOA11NM_x	BNG_E	BNG_N	LONG_	LAT	Shape_Leng	Shape__Area	Shape__Length	geometry	MSOA11NM_y	MSOA21CD	MSOA21NM	LAD22CD	LAD22NM	CTRY21CD	CTRY21NM	TCITY15CD	TCITY15NM
7196	7197	W02000419	Denbighshire 017	302419	380863	-3.46620	53.3157	9384.767259	2.716706e+06	9384.767259	POLYGON ((-3.45799 53.32180, -3.45803 53.32179...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
7197	7198	W02000420	Wrexham 020	336850	353606	-2.94407	53.0759	70438.643633	7.919196e+07	70438.643633	POLYGON ((-2.96987 53.13316, -2.96982 53.13316...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
7198	7199	W02000421	Ceredigion 011	271405	266660	-3.88655	52.2830	193683.750283	7.700700e+08	193683.750283	POLYGON ((-3.74490 52.50741, -3.74487 52.50740...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
7199	7200	W02000422	Cardiff 048	319643	174132	-3.15805	51.4604	22789.715621	8.381063e+06	22789.715621	MULTIPOLYGON (((-3.12022 51.38026, -3.12008 51...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
7200	7201	W02000423	Cardiff 049	318619	174985	-3.17298	51.4679	11723.546322	2.065266e+06	11723.546322	POLYGON ((-3.16140 51.47828, -3.16145 51.47816...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

	geometry	uID
0	POLYGON ((531541.552 181119.077, 531541.397 18...	0
1	POLYGON ((531367.177 181482.519, 531313.807 18...	1
2	POLYGON ((533328.512 181252.189, 533327.720 18...	2
3	POLYGON ((533194.003 181833.569, 533283.875 18...	3
4	POLYGON ((533408.038 180948.164, 533430.570 18...	4

	geometry	uID	area
0	POLYGON ((531541.552 181119.077, 531541.397 18...	0	681.109880
1	POLYGON ((531367.177 181482.519, 531313.807 18...	1	4153.462643
2	POLYGON ((533328.512 181252.189, 533327.720 18...	2	2509.550301
3	POLYGON ((533194.003 181833.569, 533283.875 18...	3	22909.381927
4	POLYGON ((533408.038 180948.164, 533430.570 18...	4	8015.503871

MSOA Map Read¶

Reverse Geocode City, State from Geometry¶

Buildings¶

Streets¶

Tessellation¶

From building footprint we can get morphological tessellation¶

Link Streets: Link unique IDs of streets to buildings and tessellation cells based on the nearest neighbor join.¶

Notes on Measures¶

Dimensions¶

Area¶

Missing MSOA by Polygon 462¶

Manual 2¶

End of Manual¶

Appendix¶

Area data¶

Notes on Geopandas¶

	geometry	nID
0	LINESTRING (531797.163 181408.047, 531794.157 ...	0
1	LINESTRING (531797.163 181408.047, 531802.260 ...	1
2	LINESTRING (531797.163 181408.047, 531793.703 ...	2
3	LINESTRING (531794.157 181410.161, 531797.163 ...	3
4	LINESTRING (531794.157 181410.161, 531787.903 ...	4

	uID	geometry	area
1890	1891	POLYGON ((533233.133 180472.251, 533233.911 18...	20841.969772
1515	1516	POLYGON ((533684.769 180681.247, 533685.651 18...	21799.984164
1604	1605	POLYGON ((533761.571 180792.217, 533761.112 18...	12045.263184
1685	1686	POLYGON ((533317.671 180516.803, 533317.045 18...	11396.011016
117	117	POLYGON ((533416.844 180570.669, 533416.442 18...	14819.723164

	geometry	uID	area	eri	index_right	nID
0	POLYGON ((531541.552 181119.077, 531541.397 18...	0	681.109880	1.001729	230	230
1	POLYGON ((531367.177 181482.519, 531313.807 18...	1	4153.462643	1.044225	684	684
2	POLYGON ((533328.512 181252.189, 533327.720 18...	2	2509.550301	1.125147	407	407
3	POLYGON ((533194.003 181833.569, 533283.875 18...	3	22909.381927	0.890865	1187	1187
4	POLYGON ((533408.038 180948.164, 533430.570 18...	4	8015.503871	0.849627	596	596

	bID	geometry
0	0	POLYGON ((531259.618 180770.803, 531248.786 18...
1	1	POLYGON ((531544.647 181527.606, 531544.286 18...
2	2	POLYGON ((531607.121 180421.019, 531597.024 18...
3	3	POLYGON ((532640.971 180764.898, 532641.496 18...
4	4	POLYGON ((532625.826 180813.963, 532625.718 18...

	geometry	uID
0	POLYGON ((559821.771 184975.475, 559824.415 18...	0
1	POLYGON ((559885.080 185120.185, 559894.057 18...	1
2	POLYGON ((559839.042 185058.051, 559844.832 18...	2
3	POLYGON ((559900.608 184966.984, 559931.914 18...	3
4	POLYGON ((559867.959 185052.251, 559892.879 18...	4

	geometry	uID
0	POLYGON ((545793.342 189910.798, 545765.725 18...	0
1	POLYGON ((545838.277 189934.961, 545845.723 18...	1
2	POLYGON ((545407.774 190526.928, 545405.558 19...	2
3	POLYGON ((544745.569 188909.599, 544728.235 18...	3
4	POLYGON ((544999.155 188620.739, 544982.296 18...	4

	geometry	uID
0	POLYGON ((535968.061 181595.423, 535983.519 18...	0
1	POLYGON ((536291.988 181218.724, 536299.901 18...	1
2	POLYGON ((536355.330 181271.590, 536354.095 18...	2
3	POLYGON ((536290.267 181375.995, 536290.478 18...	3
4	POLYGON ((535835.708 181399.502, 535837.147 18...	4